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ABSTRACT 

The dynamical friction force experienced by a massive gravitating body moving through a gaseous 
medium is modified by sufficiently strong large-scale magnetic fields. Using linear perturbation theory, 
we calculate the structure of the wake generated by, and the gravitational drag force on, a body 
traveling in a straight-line trajectory in a uniformly magnetized medium. The functional form of the 
drag force as a function of the Mach number (= Vq jcg , where Vq is the velocity of the body and Cg the 
sound speed) depends on the strength of the magnetic field and on the angle between the velocity of 
the perturber and the direction of the magnetic field. In particular, the peak value of the drag force is 
not near Mach number ~ 1 for a perturber moving in a sufhciently magnetized medium. As a rule of 
thumb, we may state that for supersonic motion, magnetic fields act to suppress dynamical friction; 
for subsonic motion, magnetic fields tend to enhance dynamical friction. For perturbers moving along 
the magnetic field lines, the drag force at some subsonic Mach numbers may be stronger than it is 
at supersonic velocities. We also mention the relevance of our findings to black hole coalescence in 
galactic nuclei. 

Subject headings: black hole physics — hydrodynamics — ISM: general — waves 



1. INTRODUCTION 

An object moving in a background medium induces 
a gravitational wake. The asymmetry of the mass den- 
sity distribution upstream and downstream from the per- 
turber produces a drag on the body, which is often re- 
ferred to as gravitational drag or dynamical friction (DF) 
force. A body in orbital motion may undergo a radial de- 
cay of its orbit due to the loss of angular momentum by 
the negative torque caused by DF drag. Chandrasekhar 
(1943) derived the dynamical friction on a massive par- 
ticle passing through a homogeneous and isotropic back- 
ground of light stars. His formula is applied to estimate 
the merger timescale of satellite systems or to study the 
accretion history of galaxies. Bondi & Hoyle (1944) con- 
sidered the problem of the mass accretion by a point mass 
M travelling at velocity Vb in a coUisional homogeneous 
medium of sound speed Cg in the limit where the per- 
turber moves at supersonic velocities relative to the am- 
bient gas (i.e. high Mach numbers). If the perturber is an 
accretor, streamlines with small impact parameter may 
become bound because of energy dissipation in shocks, 
and can be accreted to the perturber. Hence, the force on 
the perturber consists of two parts; the gravitational drag 
and the momentum accretion force. The latter contribu- 
tion may be decelerating or accelerating (Ruffert 1996). 
If the size of the perturber is larger than the Bondi- Hoyle 
accretion radius defined as Rbh = GM/[cl{l + M"^)] 
with A4 = Vq/cs, the density and velocity structure of 
the wake, at any Mach number, can be inferred analyti- 
cally in linear theory because the body produces a small 
perturbation in the ambient gaseous medium at any lo- 
cation. The gravitational drag is inferred as the grav- 
itational attraction between the perturber and its own 
wake (e.g., Dokuchaev 1964; Rephaeli & Salpeter 1980; 
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Just & Kegel 1990; Ostriker 1999; Kim & Kim 2007; 
Sanchez-Salcedo 2009; Namouni 2010). 

The studies of the gravitational drag in gaseous media 
have enjoyed widespread theoretical application, ranging 
from protoplanets to galaxy clusters. It seems to play 
a significant role in the growth of planetesimals (Hor- 
nung, PeUat & Barge 1985; Stewart & Wetherih 1988), 
the eccentricity excitation of planetary embryo orbits 
(Ida 1990; Namouni et al. 1996), the orbital decay of 
common-envelope binary stars (e.g., Taam & Sandquist 
2000; Nordhaus & Blackman 2006; Ricker & Taam 2008; 
Maxted et al. 2009; Stabler 2010), the evolution of the 
orbits of planets around the more massive stars (Villaver 
& Livio 2009), the evolution of low- mass condensations 
in the cores of molecular clouds (Nejad-Asghar 2010), the 
mass segregation of massive stars in young clusters em- 
bedded in dense molecular cores (Chavarrfa et al. 2010), 
the orbital de cay of kpc-sized giant clumps in galaxies a t 
high redshift (jlmmeh et al.ll2004l: iBOTarnaud et al.ll2007D . 
or the heating of intracluster gas by supersonic galaxies 
(El-Zant et al. 2004; Kim et al. 2005; Kim 2007; Conroy 
& Ostriker 2008). Special work has been devoted to un- 
derstand the role of gaseous DF in the orbital decay of 
stars and supermassive black holes as a result of hydro- 
dynamic interactions with an accretion flow in galactic 
nuclei (Narayan 2000). Mergers of supermassive black 
hole binaries may be accelerated on sub-parsec scales by 
angular momentum loss to surrounding gas (Armitage 
& Natarajan 2005). In particular, gaseous DF expedites 
the growth of SMBH by mergers in colliding galaxies (Es- 
cala et al. 2004, 2005; Dotti et al. 2006; Mayer et al. 2007; 
Tanaka & Haiman 2009; Colpi & Dotti 2011). 

Less developed is the corresponding theory of DF in 
a magnetized gaseous medium. As far as we know, the 
analytic estimate of the gravitational drag for a body 
moving on a rectilinear trajectory parallel to the uniform 
unperturbed magnetic field lines by Dokuchaev (1964) is 
the only work in this area. He concluded that the DF 
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force on a supersonic body is reduced by a factor that 
depends on the ratio between the Alfven speed and the 
sound speed. Since large-scale magnetic fields are ubiq- 
uitous in many astronomical systems such as molecular 
clouds (Tamura & Sato 1989; Goodman & Heiles f994; 
Matthews & Wilson 2002; Heiles & Crutcher 2005) or 
galactic nuclei, it is important to understand how the 
DF force is affected by the presence of ordered large- 
scale magnetic fields. In fact, young stellar systems and 
low-mass condensations orbiting in the potential of their 
birth clusters can interact with the surrounding dense 
and magnetized molecular interstellar medium during 
the dispersal of the cluster's gas. In the Galactic cen- 
ter, structures associated with ordered magnetic fields, 
called arches and threads, are detected in radio contin- 
uum maps (Yusef-Zadeh et al. 1984). The magnetic field 
configuration of the Galactic center has been viewed as 
poloidal in the diffuse, interstellar (intercloud) medium 
and approximately parallel to the Galactic plane only in 
the dense molecular clouds (Nishiyama et al. 2010). On 
the scale of 400 pc, fields of lOO^G have been reported 
(Chuss et al. 2003; Crocker et al. 2010). 

The importance of gaseous DF in the evolution and 
coalescence of a massive black hole binary is motivated 
by both observational and theoretical work that indicate 
the presence of large amounts of gas in the central re- 
gion of merging galaxies. During the merger of galaxies, 
the inflow of gas material towards the galactic center 
driven by tidal torques associated with bar instabilities 
and shocks will sweep up and amplify the magnetic field 
in the central region (Callegari et al. 2009; Guedes et 
al. 2011). Observations of gas-rich interacting galaxies 
such as the ultraluminous infrared galaxies (ULIRGs) 
show that their central regions contain massive and dense 
clouds of molecular and atomic gas (Sanders & Mirabel 
1996). ULIRGs are natural locations to expect very 
strong magnetic fields (Thompson et al. 2006; Robishaw 
et al. 2008; Thompson et al. 2009). 

In this paper we will study the DF in a gaseous medium 
on a body moving on rectilinear orbit in a homogeneous, 
uniform magnetized cloud. This is the simplest idealized 
extension of the unmagnetized case and is the first step in 
understanding the role of ordered magnetic fields. Previ- 
ous works have shown that, although the formulae of the 
gaseous drag force in a unmagnetized gas medium, were 
derived for rectilinear orbits in homogeneous and infi- 
nite media (Dokuchaev 1964; Rephaeli & Salpeter 1980; 
Just & Kegel 1990; Ostriker 1999; Sanchez-Salcedo & 
Brandenburg 1999; Kim & Kim 2009; Namouni 2010; 
Lee & Stabler 2011; Canto et al. 2011), simple 'local' 
extensions have been proven very successful in more re- 
alistic situations, e.g. smoothly decaying density back- 
grounds or when the perturber is moving on a circular 
orbit (Sanchez-Salcedo & Brandenburg 2001; Kim & Kim 
2007; Kim et al. 2008; Kim 2010). As a useful starting 
point for understanding the role of a large-scale magnetic 
field, we also consider that the unperturbed medium is 
homogeneous and uniformly magnetized. A discussion 
on the DF in other initial force-free configurations will 
be given in a separate paper. Even in the simple case 
of a uniformly magnetized medium, the magnetic field 
produces qualitatively new phenomena. 

The paper is organized as follows. In ^ we discuss the 
basic concepts on the ideal problem of a particle travel- 



ing at constant speed through a uniform gas, both in 
the purely hydrodynamic case and when the plasma is 
magnetized. In ^ we outline the linear derivation for 
calculating the steady-state density wake generated by 
an extended body moving along the magnetic fields, give 
an analytical solution of the problem and compare it with 
previous work. The time-dependent linear perturbation 
theory is presented in 21 ^ describes the structure of 
the resulting wake and evaluate the DF force as a func- 
tion of Mach number, for different angles between the 
direction of the pcrturber's velocity and the direction of 
the magnetic field. In f}6l we summarize our results and 
discuss their implications. 

2. DYNAMICAL FRICTION IN GASEOUS MEDIA: BASIC 
FORMULAE 

2.1. Unmagnetized medium 

Under assumption of a steady state, Dokuchaev (1964), 
Ruderman & Spiegel (1971) and Rephaeli & Salpeter 
(1980) derived the drag force on a point mass M mov- 
ing at velocity Vq on a straight-line trajectory through a 
uniform medium with unperturbed density pQ and sound 
speed Cs- For subsonic perturbers {A4 = Vq/cs < 1, 
where Ai is the Mach number), these authors found that 
the drag force is zero because of the front-back symme- 
try of the density distribution about the perturber. For 
the steady-state supersonic case, the drag force takes the 
form 

Fdf = 779 In A, (1) 



U2 



where A — 



max/ ' mill ; 



being rmax and rmin the maximum 
and minimum radii of the effective gravitational inter- 
action of a perturber with the gas. For extended per- 
turbers, r^in is its characteristic size, whereas for point- 
like perturbers r^iin is of the order of the Bondi-Hoyle 
radius Rbh (Canto et al. 2011), as defined in ^ 

Using a time-dependent analysis in the unmagnetized 
case, Ostriker (1999) found that (1) the force is not zero 
at the subsonic regime because, although a subsonic per- 
turber generates a density distribution with contours of 
constant density corresponding to ellipsoids, there are 
always cut-off ones within the sonic sphere that exert a 
gravitational drag, and (2) rmax increases with time in 
the supersonic case. More specifically, she found that the 
Coulomb logarithm is given by: 



lnA=-lnU-^)-^, 



for < 1 and t > r,^in/{cs — Vq), and 



(2) 



(3) 



for > 1 and t > rinin/(Vb — Cs)- The perturber is 
assumed to be formed at i = 0. The transition be- 
tween the subsonic to the supersonic regime is smooth 
without any divergence at a Mach number of unity (see 
Fig. 3 in Ostriker 1999). Sanchez-Salcedo & Branden- 
burg (1999) tested numerically that Ostriker's formula 
is very accurate for non-accreting extended perturbers. 
In many astrophysical situations, one needs to assign 
a softening radius to the gravitational potential which 
in turn determines rmin without any ambiguity. For a 
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body described with a Plummer model with core radius 
^soft ^ Rbh, Sanchez-Salcedo & Brandenburg (1999) 
found that — 2.25i?soft- 

For point-mass accretors, the friction force has been 
derived by Lee & Stahler (2011) in the subsonic regime 
and by Canto et al. (2011) in the hypersonic Umit. 

2.2. Magnetized medium 

The presence of a small-scale magnetic field tangled at 
scales below r„^in will change the speed of sound. For 
isotropic compression of a random magnetic field, the 
effective sound speed is (cj + |c^)^/^ (e.g., Zweibel 2002), 
where Ca is the Alfven speed of the random small-scale 
component of the magnetic field. Therefore, in order to 
include the effect of a small-scale magnetic field, one has 
to replace the sound speed by the effective sound speed 
in the definition of M in Eqs. ([2]) and ([3]). 

The extension of the drag force formulae is by no means 
straightforward if the gaseous medium is permeated by 
a regular magnetic field. Dokuchaev (1964) derived the 
gravitational drag force in the steady state for perturbers 
moving along the lines of the unperturbed magnetic field. 
He found that the DF drag is 
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^0 



^0 



In A, 



(4) 



Sit Vq > (c^ -I- c^)^/^, where ca is the Alfven speed of the 
regular magnetic field, and it is zero for Vq < (cl+Cj^)^^^ . 
By comparing Eqs. ([T]) and dH), we see that the drag in a 
uniform magnetized background is never larger than in 
the unmagnetized case. According to Dokuchaev (1964), 
the gravitational drag on a body with velocity Vq in a 
uniformly magnetized medium is equal to the drag on a 
body with velocity Vb/(1 — c^/Vp^)^/^ in a unmagnetized 
medium. Therefore, if one naively uses the nonmagnetic 
formulae (P)-® by replacing the sound speed Cs by the 
magnetosonic speed (c^ -I- c^)^^^ would yield to wrong 
results. In the next Section, we will show, however, that 
the paper of Dokuchaev (1964) contains an error and it 
is not true that the drag force in the magnetized medium 
case is always smaller than in the unmagnetized case. 

As Ostriker (1999) demonstrated in the field-free case, 
the steady state result found by Dokuchaev that the 
net force is zero at Vq < (c^ -I- c^)^/^, because of the 
front-back symmetry of the density perturbation in the 
medium, may be misleading. It is also unclear how In A 
varies in time for the magnetized supersonic case. More- 
over, it left unexplored how the drag force depends on 
the angle between the velocity of the perturber and the 
direction of the magnetic field. Before addressing these 
questions, however, it is still worthwhile finding analyt- 
ical solutions for the perturbed steady density and the 
resulting drag force in the simplest scenario in which the 
velocity of the perturber and the magnetic field are par- 
allel. Such a exact treatment will allow us to gain insight 
into more complicated situations. This will be done in 
the next Section. 

3. AXISYMMETRIC CASE: VELOCITY OF THE 
PERTURBER PARALLEL TO THE DIRECTION OF THE 
MAGNETIC FIELD 

We consider a gravitational perturber moving on a 
straight-line at constant velocity in a medium with un- 



perturbed density po and thermal sound speed c^. In 
the absence of magnetic fields, the linearized equations 
of motion can be reduced to a nonhomogeneous wave 
equation for the relative perturbation {p — po)/po (e.g., 
Ruderman & Spiegel 1971; Ostriker 1999). Once a uni- 
form magnetic field, Bq, parallel to the direction of per- 
turber's velocity is included, Dokuchaev (1964) showed 
that the relative perturbation obeys an equation of fourth 
order in t and solved it using a double Fourier-Hankel 
transformation. As it will become clear later, we prefer 
to describe the evolution of the system through wave- 
equations because it facilitates the physical interpreta- 
tion of the problem and because the contact with the 
analysis of Ostriker (1999) is easier. In addition, the 
extension of the equations for a case where the angle be- 
tween Bq and the velocity of the perturber is arbitrary, 
becomes straightforward in our approach. 

3.1. Perturbed density distribution 

We study first the completely steady flow created by a 
mass on a constant-speed trajectory parallel to the lines 
of the unperturbed magnetic field Bq — BqZ. To do this, 
consider a particle at the origin of our coordinate system, 
surrounded by a gas whose velocity far from the particle 
is Vq — —VqZ, with Vq > 0. We will further assume 
that the gas evolves under flux-freezing conditions. Our 
analysis begins with the linearized MHD equations to 
describe the medium's response to the perturber's pres- 
ence p = pq + p', V = Vq + v' and B — Bq + B' , in a 
stationary sate [d/dt = 0) 



Po V -v' + Vq- Vp' = 0, 



(5) 



(Vo-V)-y' 



clVp' 



- V$ 



1 



■(VxB')xBo, (6) 



Po 47rpo 
V X {VqX B') + V X {v' X Bq) = 0, (7) 

V • B' = 0, (8) 

where $ is the gravitational potential created by the per- 
turber. The Poisson equation links the potential with the 
density profile of the perturber pp: 



inGpp. 



(9) 



The Lorentz force, which provides the magnetic back- 
reaction on the fiow pattern, is given by 



{VxB')xBo = Bq 



dBL dB' 



dz 



dx 



X — 



dB' 



dy 



dz 



y 



Hence, the divergence of the Lorentz force is 
V • [(V X B') X Bq] = 

( d^B' d^B' d^B' 



(10) 



= Bn 



dxdz 



dy'^ dydz 



= -BnV^B'. 



(11) 



In the last equality we have used that V • B' = 0. By 
substituting equations ^ and (fTTj) in the divergence of 
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the equation of motioifl, we have 



Pa dz^ 



Po 



Bp 



V^B^ (12) 



By comparing the second and third terms in the right- 
hand side of the equation above, we see that, formahy, 
the magnetic back-reaction term V^B'^ is mathematically 
equivalent to having an external potential term. How- 
ever, whilst $ is known (Eq. B'^ is coupled to the 
fluid motions through the flux- freezing equation ([7]). 

Next, we need an independent equation for B'^ to close 
the system. This can be accomplished using the third 
component of the induction equation ([7]), which has the 
form 

B„fv-.'-2!i).roS. (13, 



dz 



dz 



Our strategy is to eliminate v' in Equation (|13p . The 
third component of the equation of motion can be 
written as 



dz 



Pq dz 



(9$ 

dz 



(14) 



This equation does not depend explicitly on the frozen-in 
magnetic field because the z-component of the Lorentz 
force vanishes in the linear approximation (see Eq. I10|). 
Substituting Eqs. ^ and in Eq. (HH), we obtain 
the desired equation 



Bo 



]_df/_ 

Po dz 



1 dp' 1 9$ 



poM'^ dz 



I/q^ dz 



dz 



(15) 



where we recall that M = Vo/cg is the (sonic) Mach num- 
ber. Once again, the magnetic term dB'^/dz is formally 
identical to d^/dz, but some caution should be used 
when interpreting it; the z-component of the Lorentz 
force is not dB'^/dz but zero. 

Equations (IT5|) and ([T5|) constitute a system of two 
coupled linear differential equations for p' and B'^ which 
may be solved once we have chosen suitable bound- 
ary conditions. Defining the dimensionless perturbations 
a = p' / Po and (3 = B' /Bo, the equations to solve are: 



M 



, d^a 
dz^ 



1 



{M^ - 1) 



da 1 d<P 

dz c? dz 



,di3^ 

dz 



(16) 



(17) 



where T = ca/cs and ca the Alfven speed in the un- 
perturbed medium. In the limit of vanishing magnetic 
field, = T = and Equation (|16l) reduces to that of 
the wake of a body in a unnmagnetized medium (e.g., 
Sanchez-Salcedo 2009). 

For a point-like perturber of mass M, an analytical so- 
lution can be derived for the density enhancement, veloc- 
ity and magnetic fields in the wake. In order to calculate 

^ The curl of the equation of motion provides a relationship 
between the vorticity u: and the current density J': 



Vo 



dui _ Bo dJ' 

dz cpo dz 



In linear theory, the baroclinic term vanishes and the Lorentz term 
is the only able to generate vorticity, even if the gravitational force 
is irrotational. 



the dynamical friction force exerted on the body, we only 
need the gas density enhancement in the wake, which is 
derived in Appendix and is given by 



a{R, z) = 



X{l~ri)GM 



1 



(18) 



where R — + is the cylindrical radius and 

V = {cA/Vof = {T/Mf, (19) 
£.^1 + {1-M-^)T^ = 1-T] + T\ (20) 

7 =1- — 

and 



(21) 



1 if < A^crit; 

2 if TWcrit < < min(l, T) and z/R > 

1 if min(l,T) < < max(l,T); 

2 if A^ > max(l,T) and z/i? < -|7|; 
otherwise. 



Here, the critical Mach number is defined as 



M„it = (1 + T- 



-1/2 



(22) 



(23) 



Because of the linear-theory assumption. Equation (|18p 
is properly vahd only for (z^ -|- 7^i?^)^/^ » (1 — 
ri)GM/{£^c^). The nonmagnetic steady-state solution for 
density in the wake past a gravitating body is recovered 
when T = 0. 

For clarity, it is convenient to distinguish four intervals 
depending on the value of the Mach number of the body: 
M < Alcrit (case or interval I), Alcrit < A^ < min(l, T) 
(case or interval II), min(l,T) < Ai < max(l,T) (case 
III) and Ai > max(l,T) (interval IV). 7^ is a posi- 
tive number in cases I and III, whereas it is negative 
in cases II and IV. In the latter cases where 7^ < 0, 
the density perturbation a vanishes at some spatial lo- 
cations. In case IV, for instance, a outside the cone 
defined by the surface z — — |7|i? is actually zero. Turn- 
ing to Eq. (jlSp . we see that the magnetic perturbation 
B'^ in these regions does not vanish but obeys the fol- 
lowing relation, dB'^/dz = —{Bo/VQ)d^/dz. Now, from 



Eq. (I14p . the axial component of the velocity satisfies 
a similar equation, dv'^/dz = {l/Vo)d^/dz. Using the 
fact that V • t>' = in regions of constant density (see 
Eq.[S]), it is simple to show that B' is parallel to v' in re- 
gions where a = 0, and thus the magnetic configuration 
is force-free in these zones. 

Once p/po is known, the gravitational drag can be com- 
puted; this will be done in Section Nevertheless, in 
order to gain more insight into the physics of the wake, we 
will describe the morphology and structure of the steady- 
state wake in the next Section. 

3.2. Physical interpretation 

Consider first subsonic perturbers. In the limit Ai — !• 
0, we have 7^1,$^ -T^/M^, (1 - r;) ^ ~T^/M^ 
and A = 1. Therefore, the density enhancement is 
GM/{c'^r), which is T-independent, and corresponds 
to the linearized solution of the hydrostatic envelope, 
p/po = exp [GM/{clr)'\, around a stationary perturber 
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(e.g., Ostriker 1999). In this case, the magnetic field hnes 
remain straight and the whole magnetic configuration is 
force-free. 

The surfaces of constant density for subsonic per- 
turbers may be either ellipsoids or hyperbolae, depend- 
ing on the Mach number. At < AicHt it holds that 
7^ > 1 and the isodensity surfaces are ellipsoids elon- 
gated along the trajectory of the perturber with eccen- 
tricity e — Ai/ V^ICT- This is in sharp contrast to what 
happens without any magnetic field where the ellipses are 
elongated along R for perturbers at any subsonic Mach 
number (7^ < 1). Therefore, the existence of ellipsoidal 
density isocontours elongated along z is a clear signature 
of curved magnetic fields. 

In the nonmagnetic subsonic case, < at z > and 
v'^ > at any z , signifying that the incoming fiuid is veer- 
ing towards the perturber, but then turning away again 
once it passes the body. The result oiv'^>0 implies that 
the gas is being dragged by the gravitational pushing of 
the body. In order to understand the morphology of the 
wake in a magnetized medium, in the following we cal- 
culate v'j^ and v'^ for a perturber with Mach number in 
the interval I. 

From Eq. ([H]) and using the result for a in Eq. ([T5)) . 
we find that 



at not extremely large z distances, is 



Vo 



dv' 



Ml - 11) 



GMz 



dz 

which leads to 
GM 



GMz 



£, (z^ + R^-f^f/^ (z2 + i?2)3/ 



A(l-,;) 



Vq V(^^+i?2)l/2 ^(^2,^^2^)1/2 



(24) 



(25) 



It is simple to show that > in case I, regardless the 
value of T. 

The radial component of the velocity can be found us- 
ing Eq. ©, 



1 dRv'j^ 
R dR 



dv' da 
dz dz 



7W2 



{l-M')— + 
dz 



da 1 a$ 



dz 
(26) 

Since we already know a and $, this equation can be 
solved to obtain the radial velocity: 



GMz 



1 



A(l-7W2)(l-r?) 

VoR \{z^ + R^)^/^ -f^^iz^ + R^-f^y/^ 

'(27) 

From the equation above, it follows that, in case I, > 
at the head of the perturber, regardless the magnetic field 
strength. Thus, a parcel of fluid in the upstream region 
circulates around the perturber, reaching its maximum 
i?-value at z = and then turning back again. Since 
frozen-in magnetic field lines are dragged by the gas, the 
upstream magnetic field lines are decompressed radially, 
resulting in arched magnetic field lines with negative B'^- 
values. In fact, integration of Eq. ([15]), gives 



B' = 



GMBo 

T/2 



1 



{z^ + R^y/-^ ^{z'^ + R^-f^y/^ 



(28) 

This clearly states that B'^ < (here we assume Bq > 0). 
Thus, an anticorrelation between density and B'^ arises. 

At Mach numbers close to A^crit, that is, Ai = A^crit — 
e with e a very small positive number, the density profile. 



a{R, z) 



T^'^GM/cl 



2e{l + T^y/^R' 



(29) 



Hence, the density enhancement is large and its z- 
gradient very small. 

At subsonic Mach numbers in the interval Atcrit < 
M < min(l, T), the surfaces of constant density exhibit 
no front-back symmetry. The isodensity contours corre- 
spond to hyperbolae in the z — R plane, as occurs for 
supersonic perturbers in the absence of magnetic fields, 
but now the density perturbation is null in the rear Mach 
cone and is non-vanishing in a modified Mach cone lead- 
ing the perturber. The physical reason is as follows. The 
density perturbation a is non-positive at any location 
because the streamlines diverge at the front cone (?? > 1 
and ^ > 0). a at the edges of the cone is minus infinity 
for a point mass. In the front cone, dv'^/dz > 0, mean- 
ing that the fiow in that region is being accelerated by 
the inward net pressure force. Across the edge of the 
modified Mach cone, there is a rapid rise in pressure and 
density, and the gas velocity quickly slows. In fact, the 
causality criterion used in Appendix lAl is tantamount to 
selecting the solution in which a rapid fiow is slowed in a 
short distance, as occurs in shock waves. Indeed, we will 
find numerically further below that the system adopts 
this solution (i i5.1|) . 

In case III, the body moves at intermediate velocities, 
i.e. either in the range Cs < Vq < ca (if T > 1) or in the 
range ca < Vq < Cg (if T < 1). It is only in this case that 
7 lies in the range < 7 < 1 and, therefore, the ellipsoids 
are flattened along z. Remind that in the unmagnetized 
background, a subsonic perturber also produces ellipti- 
cal density distributions flattened along z (e.g., Ostriker 
1999); the latter is a particular case oi ca < < Cg- 
For Cs <Vq < Ca (which requires that T > 1), however, 
the density perturbation is negative {a < 0). In this 
case, the thermal pressure decreases towards the body. 
Therefore, ahead of the perturber the gas is accelerated 
along the z-direction by the gravitational force plus the 
pressure gradient (it holds that dv'^/dz > at z > 0). A 
negative a is a consequence of the action of the pressure 
gradient plus the reduction of the radial convergence of 
the flow due to the presence of the ordered magnetic field 
that preferentially allows motions along z, which lead to 
an accelerated flow falling towards the perturber. In the 
radial direction, the magnetic fleld lines are compressed 
at z > because vr < upstream of the perturber. 

From Eqs. p^ - (|22p we see that if the motion of the 
perturber is both supersonic and super- Alfvenic, which 
corresponds to case IV, the density disturbance is con- 
fined to a rear cone defined by the condition z < — |7|i?. 
In this regime, the surfaces of constant density within the 
wake correspond to similar hyperbolae in the z — i? plane, 
at the rear of the body, with eccentricity e = Ai/^/^. 
From Eq. 1^^, it is simple to show that ^ > 1 in this case. 
Therefore, given a sonic Mach number A4, the angular 
aperture of the cone is larger in the presence of mag- 
netic fields. In analogy to the unmagnetized medium, 
one could define the effective speed of propagation of 
the disturbance as Vp = y/^Cs in case IV. At large Mach 
numbers (say A4 3> \/^), the cone is very narrow (e ^ 1) 
and, as expected, the propagation speed coincides with 



6 



SANCHEZ-SALCEDO 



the magnetoacoustic velocity, Vp ~ -y/cf + cj. Hence, 
at these large Ai values, the stationary flow is similar 
to that in the nonmagnetic case but replacing the sound 
speed by the magnetosonic speed. 

Now consider case IV but when the perturber moves 
at the same velocity as the effective speed of the distur- 
bance, so that Vq = Vp 01, equivalently, M = In the 
non-magnetic case, this condition corresponds to ^ = 1 
and, therefore, the velocity of the body is in resonance 
with the sound speed in the medium. One could naively 
think that, at A4 = y^, the response of the medium 
is maximum because of the resonance between Vq and 
Vp. This is not true for T > 1 because Ai = y/£^ implies 
M^T (using Eq.[20l), = 1 (from Eq.^ and, thereby, 
a = 0. We learn that a mass moving at the Alfven speed 
in a medium with T > 1, does not generate any density 
disturbance in the ambient gas because the velocity field 
of the stationary flow is divergence- free (V ■ v' = 0). 

3.3. A comparison with Dokuchaev (1964) 

As already mentioned, Dokuchaev (1964) calculated, 
for the first time, the properties of the wake created by 
a star moving along the field lines, by treating it as a 
linear perturbation. His analysis started from the time- 
dependent linearized equations of magnetohydrodynam- 
ics, including a source term Q in the continuity equa- 
tion, representing the gas replenishment by the star. Al- 
though he used the time-dependent equations, he tac- 
itly assumed that the object's gravitational field is active 
since t = —oo, so that the wake is in a steady state. For 
the case without mass injected by the star ((3 = 0), the 
physical stand points used by Dokuc haev (1964) are ex- 
actly the same as those adopted in except that he 
chose a reference frame in which the unperturbed back- 
ground gas is at rest. 

Dokuchaev (1964) found closed fourth-order differen- 
tial equations for p and the radial component Vr of the ve- 
locity. Through Fourier-Hankel transformations, he was 
able to solve the differential equation for p. He found a 
similar expression for p as that given in Eq. (jl8p but he 
failed to separate correctly the different intervals for A 
(Eq. [22]) and the intervals at which the isocontours are 
ellipsoids or hyperboloids. In particular, he claimed that 
the isocontours are ellipsoids at any Mach number below 
(c^ + c\y/'^/cs = (1 + T^)^/^, which is misleading. 

3.4. Gravitational drag force in the axisymmetric case 

Once we have the gas density enhancement a{r) in the 
ambient medium, we can calculate the gravitational force 
exerted on the body by its own wake: 



where 



^ DF = 



dz dR Ra(r) 



(30) 

The net drag is zero when the isodensity contours are 
ellipsoids, i.e. when 7^ > 0, because the wake exhibits 
front-back symmetry. At values 7^ < 0, however, the 
region of perturbed density is confined to a cone and the 
drag force is nonvanishing. Evaluating the integrals in 
spherical coordinates {R = r sin and z = r cos 6) , and 
using the variable p defined as p = cos 9, the drag force 
can be expressed as: 



DF 



(1 - V) — — — ^ 



1/2 



dr , 
—z, 

r 



dp 



(32) 



As already said, the drag force is nonzero in cases II and 
IV, where 7^ < and thus ^ > 0. In case II, A = 2 
for all p between /xiower = (A^^ — iY^"^ /M and 1, so 
that 1=1. In case IV, A = 2 for all p between — 1 
and Supper = —{M^ — S)^^"^ /M and thus I = —1. Since 
1 — 77 < in case II, we can write (1 — 77)1 = |1 — 77I and 
the resultant expression for the force is: 

( -\l-T]\T\nAz ff7W>max(l,T) 
Fdf = < or Xcrit <M< min(l, T); 

[ otherwise 

(33) 



where 



-F = 



AnG^M^pn 



(34) 



(31) 



and In A is the Coulomb logarithm. Note that the den- 
sity diverges in the wake at Mach numbers close to A^crit 
because ^ — 0. However, the drag force is finite because 
the opening angle of the cone becomes very narrow. Still, 
the drag force peaks at Mach numbers close to A^crit be- 
cause the factor |1 — ri\/VQ in the formula for a increases 
when A4 decreases. 

Figure [1] shows the DF force felt by the body at 
t = lOOrniin/cs, as a function of the Mach number and 
for different values of T . In analogy to the unmagnetized 
case, we take A = Vb^/^min = lOOVo/cj. Dokuchaev 
(1964) claimed that the drag force is nonzero only at 
M > {1 + T2)1/2 (see This is incorrect. For in- 

stance, for T = 1, the drag force is different from zero at 
Ai > A^crit = 0.70. If the ratio between the Alfven and 
sound speeds is of = 2, the DF force is nonvanishing 
in the intervals 0.816 < A^ < 1, and M > 1.41. In fact, 
there exists always a subsonic velocity range at which the 
drag force is nonzero. 

As long as T ^ 0, the DF force has two local max- 
ima; one located at A/fcrit and the other one at Ai = 
max(l, a/2T). The drag force strength at A^crit increases 
with T, whereas the drag force at the second local maxi- 
mum decreases with T. As Figure [T] clearly shows, at low 
T- values, the width of the interval with Fjyp ^ around 
A^crit becomes very narrow. For instance, the width of 
that interval is only of 4 x 10"'^ for T = 0.2. Hence, 
the drag force at subsonic values is irrelevant for astro- 
physical purposes when the Alfven speed is sufficiently 
small as compared to the sound speed. For T > 0.4, 
the drag force at the local maximum A^crit is always 
larger than the drag force at the other local maximum 
A4 = max(l, \/2T). In the particular case of T = 1.41, 
the drag force is a factor > 6 stronger in the interval 
0.816 < A^ < 1 than it is at A^ > 1. 

4. TIME-DEPENDENT EQUATIONS 

The steady state analysis in the axisymmetric case pre- 
dicts zero drag force at certain Mach numbers because 
the perturber is surrounded by complete ellipsoids that 
exert no net force. As Ostriker (1999) demonstrated in 
the field-free case, the time-dependent analysis in which 
the body is dropped suddently at t = allows to capture 
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Fig. 1. — Gravitational drag force as a function of Mach number, at t 
axisymmetric case, for different values of T. 



lOOrmin/cs, as predicted by the steady-state linear-theory in the 
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Fig. 2. — Color map of the density p/po (left panel) and magnetic 
field B'^/ B(j (right panel), in the (i?, z)-plane for a case with Vq = 
and T = 1.41, in a natural logarithmic scale. 



the asymmetric density shells in the far field which exert 
a gravitational drag on the body. Other advantage of 
the time-dependent approach is that, contrary to what 
happens when assuming steady-state, the ambiguity in 
the definition of the maximum cut-off distance r^ax is 
fixed. 

Without loss of generality, it is convenient to use the 
gas frame of reference in which the ambient gas is initially 
at rest, the initial magnetic field is along the z-axis and 
the body moves with velocity Vb^^y -I- Vq^zZ. The first 
order continuity equation is 



(35) 



the MHD Euler equation 

dv' _ clVp' 
'dt ~~ 



- V<I> 



1 



Pa 47r/9o 
and the induction equation: 
dB' 



(V X B') X Bo, (36) 



dt 



= V X (t)' X Bo)- 



(37) 



The medium initially uniform will respond to the grav- 
itational pull of the body through the emission of fast 
and slow Alfven waves and sound waves. In the follow- 
ing we will manipulate the above equations to obtain a 
closed system of two differential equations for p' and B'^ 
in analogy to the steady-state case. 
Using Eq. (fTTI) in the divergence of Equation (p6| 

Po 



dt 



Bo 



-V^B'. 



By substituting Ec^. (1551) into Eq. 

V2$ + 



1 d'^p' cl 



Po ot^ Po 
In terms of a and ji^, it yields 



47rpo 

we obtain 

Bo 



(38) 



47r/9o 



a a 



(39) 



(40) 



Here, the magnetic effect on the density perturbation 
appears as a inhomogeneous term. We may recover the 
classical non- magnetic equation for a by taking ca = 0- 
On the other hand, the third component of the induc- 
tion equation (Eq. I37p implies: 



dBi 
dt 



- -Bn 



dx 



dy 



(41) 



dz 



(42) 



Equations (|3S]) and (gT]) give 

dt \po dt 

From the third component of the equation of motion 

dv^^_^dp^_d^ (43) 
dt po dz dz ' 
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Fig. 3. — Color map of the density p/poi in the (/?, 2)-plane, for the cyhndrical case with M = 0.75 (upper panels) and for M = 0.9 
(lower panels), in a natural logarithmie scale. Both mo dels have T = 1.41, implying A^crit = 0.816. Therefore, M = 0.75 falls into the 
interval I, while M = 0.9 lies in the interval II (see ^3.111 . To easy comparison, the density map in the steady-state for a perturber seated 
at i? = z = is shown in the left panels. The central and right panels display the density in the wake at two snapshots, when the perturber 
is dropped suddenly at t = at the origin of the coordinate system. The time of the snapshots is given in the lower right hand corner of 
each panel in units of tcross ■ 



we know that 



d_ 

di 



dz 



dz 



(44) 



(00 dz"^ 

Inserting Eq. (UH) into the temporal derivative of 
Eq. (gSI), one finds 



1 d'^B', I 9V' 



In dimensionless form: 



92$ 



(45) 



Putting together, the equations pOj) and (l46l) to solve 
can be written as 



(47) 

(a + $), (48) 



where $ is the gravitational potential in units of 

(i.e. = '^/c^s) and we have used the Lorentz invariant 
D'Alembertian □, defined as: 



1 w2 

V 



(49) 
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The first equation (Eq. |47|) governs the evolution of the 
density in the presence of a gravitational potential and 
magnetic fields. The second equation (Eq|35]) describes 
the evolution of a frozen-in magnetic field when the gas 
is subject to pressure gradients and to an external gravi- 
tational potential. The inhomogeneous term in Eq. (j48p 
does not have ^-derivatives because gas motions in that 
direction does not compress, stir or stretch the back- 
ground magnetic field. The resulting equations (j47|) and 
conform to a set of two coupled non-homogeneous 
wave equations. For a point-mass perturber, it is sim- 
ple to find a in the Fourier-Laplace space, di{k^uj)^ but 
the inverse Fourier-Laplace integral cannot be given in a 
closed analytic form. 

In order to gain physical insight, consider first a two- 
dimensional example. If $ = $(a;,y), that is, if the per- 



turber is an infinite cylinder with a certain radial den- 
sity profile pp — Pp{R), then j3z = a because of the flux- 
freezing condition, and a satisfies a simple wave equation 
with magnetoacoustic speed: 



92 



1 



(50) 



The physical reason is that motions are always perpen- 
dicular to the frozen-in magnetic field lines. Magneto- 
hydrodynamical equilibrium is reached within the mag- 
netosonic cylinder. At a later stage, Parker instabilities 
can develop (Sanchez-Salcedo & Santillan 2011). 

In the purely hydrodynamical problem, the equation 
governing the evolution of a is Uga = V^*!". If the per- 
turber is a point source, we have U^a = {A-kGM / cl)5[x— 
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Fig. 5. — Same as Fig.[3]but for M = 1.7. Again T = 1.41. This Mach number lies in the interval IV. 
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Fig. 6. — Distributions of the perturbed density a (solid lines) and the ^-component of the perturbed magnetic field (dashed lines) along 
a cut at i? = at two different times, for the same model as that shown in Fig.[5](A^ = 1.7 and T = 1.41). 
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Fig. 7. — Temporal evolution of the gravitational drag force in the axisymmetric model, at six different Mach numbers between 0.55 and 
2, for T = 1 (left panel), and for T = 1.41 (right panel). The number on eax;h curve is the Ma<;h number M. 
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Vo,xt)S{y— Vo,yt)S{z — Vo^zt)H{t). Hence, the density re- 
mains unperturbed outside the causal region for sound 
waves (see Ostriker 1999). In a magnetized medium, 
however, the situation is different because Equation (|48p 
for the perturbed magnetic field (5^ has a source term 
[d^^/dx^ + d^^/dy^) which does not vanish even out- 
side the causal region for magnetosonic waves. 

In the next Section we will solve the coupled wave- 
equations numerically. To do so, the perturber gravita- 
tional potential will be modeled by a smooth core Plum- 
mer potential: 

^, ^ GMnjt) 

V^' + iy- Vo,ytr + - Vo,.ty + b:^,, ' 

.(51) 

where -Rsoft is the softening radius and H is a Heaviside 
step function. Stellar and globular clusters can be ac- 
curately described by Plummer potentials. These type 
of models were also used in Sanchez-Salcedo & Branden- 
burg (1999, 2001), Kim & Kim (2009), and Kim (2010) 
to study DF. 

5. RESULTS 

The coupled inhomogeneous wave equations (|47p and 
(j48|) were solved using a finite difference scheme in a uni- 
form grid. The scheme is second order in space and third 
order in time. The temporal algorithm was described 
in Sanchez-Salcedo & Brandenburg (2001). Calculations 
start with a uniform background density and magnetic 
field, and the body is initially placed at the origin of the 
coordinate system with velocity Vb,j,y -I- Vq^zZ. For the 
axisymmetric case, which occurs when Vo^y = 0, the cal- 
culations were carried out on a two-dimensional (i?, z)- 
plane in cylindrical symmetry. In the general three- 
dimensional case {Vo^y ^ 0), the variables a and /?z are 
symmetric about the plane a: = 0. Hence, we considered 
a finite domain with x € [0,Lj^^niax] and used symmet- 
ric boundary conditions at x = 0, and outflow boundary 
conditions in the other five caps of the computational do- 
main. However, the size of the domain was taken large 
enough to ensure that the perturbed density and mag- 
netic field do not reach the boundaries. 

As a test of the algorithm, we studied the convergence 
of homogeneous wave modes by perturbing a uniform 
background medium. We further tested convergence of 
our models for several resolutions and found that four 
zones per i?soft suffice to have converged results. 

We take i?soft, c^, and tcross — -Rsoft /cs as the units of 
length, velocity, and time, respectively. A model can 
thus be specified with four dimensionless parameters: 
GM/{cj, Rsoft), M, T and e = atan(K),a/Vb,z). 6 is the 
angle between Vq and Bq. Fixed M, 6 and T, the vari- 
ables a and /?z depend linearly on GM/{clRsoh)- Hence, 
in our calculations we always take GM/{c'^Rso{t) = 
0.01/3 and explore how the density and the magnetic 
field in the wake depend on the other three parameters. 

5.1. Axisymmetric case 

We first run models with the magnetic field terms swich 
off and compare the density enhancement and the gravi- 
tational drag with previous linear calculations in Ostriker 
(1999) and Sanchez-Salcedo & Brandenburg (1999). We 
found excellent agreement, backing up our numerical 



model. In the following, we will present results for a body 
moving along the field lines of the unperturbed magnetic 
field, which corresponds to = 0. 

The simplest scenario occurs when the gravitational 
perturber is dropped at < = at rest {Vq — 0). As dis- 
cussed at the begining of ^'S.2l the steady-state density 
distribution is identical as that without any magnetic 
field. However, the inital relaxation stage and the far 
density distribution are sensitive to magnetic effects. In 
fact, while the problem has spherical symmetry at any 
time in the purely hydrodynamical case, this symmetry 
is broken in a magnetized medium because the magnetic 
field dictates a preferential direction. Figure [5] shows 
maps of density and B'^ for a case with T = 1.41. We 
see that the density distribution in the vicinity of the 
body is indeed spherically symmetric and the magnetic 
field takes essentially its initial value, implying that this 
part has reached hydrostatic equilibrium. However, there 
are two symmetric underdense regions along the z-axis in 
the outer parts. Physically, the origin of them is that the 
magnetic field reduces the flow convergence toward the 
symmetry axis because magnetic forces mainly affect the 
radial component of the velocity vjf . This loss of radial 
convergence produces a wave with negative density en- 
hancement (a < 0) but a positive magnetic enhancement 
{(3z > 0) because of the compression of the magnetic field 
lines in the radial direction. 

Figure [3] shows snapshots of the density at the (i?, z) 
plane for T = 1.41 and two subsonic velocities {M = 0.75 
and Ai = 0.9). In the time-dependent analysis, both 
cases present a region of negative density enhancement 
(i.e. a < 0) at the head of the perturber. As Ai increases 
from to 0.9, the underdense region at the head of the 
body remains and gets deeper, while the underdensity in 
the downstream region becomes less pronounced. This 
is simply consequence of the Doppler effect; gradients 
become steeper upstream (remind that this also occurs 
in the purely hydrodynamic case). By comparing the 
density at two times (the central and the right panels of 
Figure [3]), we see that the evolution of the density looks 
self-similar. 

At M = 0.75, the steady-state analysis predicts a null 
drag force because of the front-back symmetry in the 
wake (see upper left panel of Fig. [3]). In the finite-time 
case, complete ellipsoids are visible only in the vicinity of 
the body. For instance, a.t t = 173.6icross (with T = 1.41 
and Ai = 0.75), the isodensity contours are not longer 
ellipsoids at distances beyond ~ 8-Rsoft from the body's 
center. For comparison, in the absence of magnetic fields 
(T — 0) and at ^Vf = 0.75, ellipsoids within a radius of 
43i?soft are complete a.t t = 173.6icioss- This difference 
is a consequence of the coupling between a and f3z ■ At 
z = 136i?soft and R — (i.e. in the symmetry axis), 
a steep front separates the fluid into two regions; one 
with a > and /3z < 0, from another with a < and 
Pz > This front leads to a deceleration of the gas, which 
expands radially. 

At Ai = 0.9, a cone of negative density enhancement 
is located at the head of the perturber (see Fig. as 
that predicted in H3.ll and a region of positive density 
enhancement at the rear. Because of the back-reaction 
of the magnetic field, the isodensity contours of the tail 
are not incomplete ellipsoids at all. Note that the body 
is at the apex of the cone, whereas the overdensity is 
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Fig. 8. — Gravitational drag force for the time-dependent axisymmetric models against the Mach number at t = 40tcross (left panel) and 
at t = 200tcross (right panel). Symbols correspond to the numerical models. The solid curves plot Ostriker's formula, while the remainder 
curves draw the drag force as predicted by Equation I I33I I . adopting T'min ~ 2.25-Rsoft and Tmax — 



detached from the perturber. It is important to remark 
that, according to the analysis in i i3.4l and Figure [l] the 
maximum drag force for T = 1.41 occurs at a Mach 
number close to 0.9. 

Figure m displays density maps for T = 1.41 and two 
supersonic Mach numbers: M = 1.2 (sub-Alfvenic per- 
turber) and = 1.4 (trans- Alfvenic perturber). In the 
first case, the steady-state analysis predicts ellipsoidal 
isocontours. In the time-dependent case, however, ellip- 



soids are incomplete far enough away upstream from the 
perturber and, again, an overdensity wave at the rear 
of the perturber moves away from the body. When the 
pertuber moves at = 1.4, a teneous ellipsoidal under- 
dense envelop is still visible. 

The overdensity behind the body has the shape of a 
tulip. This tulip-shaped overdensity also appears at the 
rear of the body at M — 1.4 and at = 1.7 (see Figs. [4] 
and [5|) . A feature of the tulip-shaped overdensity wave 
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Fig. 9. — Snapshots of the density map along a cut-ofT through the {y, 2)-plane at x = 0, in natural logarithmic scale, for = 45° (upper 
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panel. The perturber was dropped at t = and moves on a rectilinear orbit in the (y, 2:)-plane. In all panels, T = 1.41. 
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Fig. 10. — Gravitational DF force as a function of Mach number 
for different T- values at t = 404cross- The perturber moves per- 
pendicular to the magnetic field lines (i.e. = 90°). The solid 
curve plots Ostriker's formula, which was derived for unmagne- 
tized media, with rmin = 2.25iJsoft a^nd ''max = VqI. To make the 
plot readable, the symbols at JV[ = 0.3, 0.55, 0.75 and 2 have been 
slightly shifted in the horizontal direction. 
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Fig . 11. — Temporal evolution of the gravitational drag force for 
T = 1. The perturber moves in a direction perpendicular to the 
magnetic field lines. The numbers given at each curve represents 
the Mach number. 

is that it has a negative magnetic enhancement. The 
tulip-shaped overdensity is a consequence of our initial 
conditions and, as expected, it is detached from the body. 
In order to illustrate the birth of the tulip-shaped wave, 
Figure [6] shows the density and magnetic perturbations 
along the symmetry axis for M = 1.7, at two early times. 
Initially, j3z increases due to the compression of mag- 
netic field lines (see the panel at t — 2.83icioss)- At the 
far edge of the tail, z ~ — 3i?soft, an underdense region 
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with positive Pz appears. Later on (see the profiles at 
t = 6.47icross), the overdensity loses gravitational sup- 
port and expands behind the body, decreasing the mag- 
netic field strength, until the magnetic pressure plus the 
magnetic tension provides sufRcient radial confinement 
to the tulip-shaped structure, allowing it to remain over 
long times. 

In Figure [5] it is simple to identify the modified Mach 
cone dragged by its point by the perturber. ki t — 
64.7fcross, the Mach cone at the rear of the body is well- 
defined. Clearly, the timescale for the development of the 
Mach cone at the rear for M = 1.7 is shorter than the 
timescale to form the upstream Mach front at = 0.9. 
In fact, Figure [3] shows that the cone for — 0.9 is not 
well developed at t = 57.7icross- 

Figure [7] shows the gravitational DF drag as a func- 
tion of time for M = 0.55,0.75,0.9, 1.1, 1.4 and 2. For 
T = 1 and M = 0.55, the drag force clearly saturates 
in ~ 30fcross- In the remainder cases the drag force in- 
creases with time. However, the drag force on a per- 
turber moving in a medium with T = 1.41, saturates for 
M = 0.55,0.75,1.1, and 1.4. This means that the DF 
force saturates to a constant value at those Mach num- 
bers that the steady-state analysis predicts a null drag 
force. 

In Fig. [S] we plot the drag force at < = 40tcross and 
t = 200icross for different values of T, together with the 
predicted force with rmin = 2.25i?soft- We see that for 
T = 0, there is a perfect agreement between the Os- 
triker formula and the inferred values, confirming the 
result that rmin ~ 2.25i?soft reported in Sanchez-Salcedo 
& Brandenburg (1999). The drag force formula given in 
Eq. (p4|) . with rmin ~ 2.25i?soft, overestimates the drag 
force in a neighbourhood of A^crit, where the drag force 
as a function of M , becomes very cuspy. As expected, the 
steady-state formula is more accurate at long timescales, 
except when it predicts a zero net force. Roughly speak- 
ing, we may say that, for the axisymmetric case, the grav- 
itational drag in a magnetized medium is always smaller 
or equal as the drag force in the unmagnetized case for 
supersonic perturbers {Ai > 1), whereas the drag is al- 
ways larger or equal in a magnetized medium at subsonic 
perturbers {M < 1). 

5.2. Magnetic field perpendicular to the direction of 
motion of the perturber 

In §5.11 we have focused on the case where the angle 
formed between the velocity of the perturber and the 
magnetic field, G, was equal to 0. In such a situation, 
the problem has axial symmetry and it is possible to find 
the analytical solution in the steady state. However, it 
is by no means clear how the gravitational drag force 
depends on 0. A visual comparison of the density wake 
structures for 8 = 0, 8 = tt/A and 8 = 7r/2 in Figures 
[3] and [5] would lead us to think that the resulting wake 
at 8 = 7r/4 is more similar to the case with 8 = 7r/2 
than to 8 = 0. In particular, we would like to stress 
the remarkably different structure of the wake for 8 = 
(Fig. El) and 8 = 7r/2 (Fig. ^ at M = 0.9. 

In this section we will discuss in detail the extreme 
case where the perturber moves perpendicular to the 
field lines. In Appendix [Bl the perturbed density is 
given in Fourier space. At subsonic Mach numbers and 
8 = 7r/2, the perturber is surrounded by a ellipsoid-like 



envelope but also presents a tail with positive and neg- 
ative a- values separated by a sharp front (see Fig. [9] for 
T = 1.41). Approaching to the upstream axis of motion 
of the body, the y-axis, the plowing up of field lines in- 
creases the total pressure. At low Mach numbers, say 
A4 = 0.3, underdense regions are now formed in the di- 
rection of the ambient field lines, which are lagged behind 
the body (remind that regions with negative a appear 
along the field lines; see Fig. [2]). At A( = 0.9, even if 
the motion is subsonic and sub-Alfvenic, a magnetic bow 
wave with sharp edges and opening angle arctan[cyi/Vo] 
is apparent in Fig. [51 Note that when we say that it is 
sub-Alfvenic, we only mean ca/Vq < 1. However, some 
caution should be used when interpreting this ratio be- 
cause the velocities are oriented in different directions. 
Since the velocity of the perturber is always orthogonal 
to the ambient direction of propagation of Alfven waves, 
the Alven speed in the direction of motion of the per- 
turber is zero and, thus, the body is always infinitely 
super- Alfvenic in the direction of motion. The morphol- 
ogy of the wake is the result of a competition between 
the gravitational focusing of gas by the perturber and 
the drainage of gas along magnetic field lines. We should 
warn here that the wake is not axisymmetric and thus 
the density map in the {x, 2/)-plane is different than the 
map in the (y, z)-plane. 

When the perturber travels faster than the magne- 
tosonic velocity Cs(l-I-T^)^/^, a magnetosonic Mach cone 
is formed at the rear of the pertuber; the entire perturbed 
density distribution lags the perturber. In the {y, z)- 
plane, the perturber creates two magnetic bow waves; 
the Alfvenic wave with opening angle arctan[cyi/Vb] and 
the magnetosonic wave with opening angle arctan[(c^ + 

The gravitational drag force is the result of the con- 
tribution of all the parcels in the domain and it is not 
possible to estimate its value just by comparing the den- 
sity structure by eye. Figure [10] shows the gravitational 
drag as a function of perturber's Mach number for dif- 
ferent values of T, together with the gravitational drag 
in the unmagnetized case using Ostriker's formula. All 
the points at Mach number larger than 0.7 lie on or be- 
low Ostriker's curve, implying that at Ai > 0.7 the drag 
force is equal or smaller than it is in the unmagnetized 
medium. At T < 0.5, the effect of including the magnetic 
field on the drag force is rather small. Interestingly, at 
T > 0.5, the strength of the drag for Mach numbers 
> (1.7 + T'^y^^ is identical as it is in the unmagnetized 
case. 

For T = 1.41, the drag force shows a plateau between 
M — 0.6 and 1.4 (see Figure [T0|). In general, the drag 
force is remarkably supressed at Mach numbers around 
^ 1 in the magnetized case as compared to the unmag- 
netized case, as long as T > 0.5 (see, for instance the 
drag at 0.9 < M < 1.4 for T = 1). In fact, the tem- 
poral evolution of the drag force is given in Figure [TT] 
for T = 1. The DF force on perturbers moving at Mach 
numbers 0.75 < Ai < 1.2 saturates asymptotically to a 
constant value. Hence, at 0.75 < < 1, the drag forces 
reach a steady-state value either the medium is magne- 
tized or not. However, we know that the unmagnetized 
drag force increases logarithmically in time for supersonic 
perturbers. This implies that the drag force may be sup- 
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Fig. 12. — Components of the gravitational DF force parallel (left panel) and perpendicular (right panel) to the direction Vo (at 
t = 40tcross) versus M, for three different values of T (T = 0.71, asterisks; T = 1, crosses; Y = 1.41, diamonds). The angle between the 
pcrturbcr velocity and the magnetic field is 45°. The key to symbols is the same as in Figure [TOl 
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Fig. 13. — Component of the DF force perpendicular to the di- 
rection of motion as a function of time for T = 1.41 and = 45°. 
The number on each curve is the Mach number. 

pressed by one order of magnitude at 1 < A4 < 1.2 in the 
magnetized as compared to the mimagnetized medium 
because the magnetized drag force saturates to a con- 
stant value. On the other hand, at < 0.75 (again 
T = 1) there is no indication that the drag force satu- 
rates, at least up to i = 40icross, implying that the DF 
force in a magnetized medium may be larger by a factor 
of a few than the drag in the unmagnetized case. 

In summary, when the magnetic fields are relevant, 
that is for T > 0.5, we distinguish three ranges. At 
high Mach numbers [i.e. M > (1.7 + T^)^/^], the drag is 
the same as in the unmagnetized case. At intermediate 
Mach numbers (1 < M < TM^.J^), the drag is highly sup- 
pressed. Finally, at low Mach numbers, the drag force is 
stronger in the MHD case than in a purely hydrodynam- 
ical medium. 

5.3. Intermediate angle between perturher's velocity and 
magnetic field. Q — 45° 

The upper panels of Figure |9] exhibit the complex mor- 
phology of the wake when 9 = 7r/4. As it is obvious from 



these panels, the gravitational drag will have two compo- 
nents: one parallel to Vq, which produces the drag and 
loss of kinetic energy by the perturber, and one compo- 
nent perpendicular to Vq, which would change the di- 
rection of Vo . Given the symmetry of the problem, both 
components lie in the [y, 2;)-plane. We will start our dis- 
cussion by considering the gravitational drag force. 

For those simulations presented in Figure [T2] with 8 — 
7r/4, the drag force saturates within t = 40icross only for 
the model with = 0.9 and T = 0.71. The maximum 
of the drag force occurs at Mach numbers near ~ M^ctAv 
At intermediate Mach numbers, say at 1.2-1.6 for T = 
1.41, the drag force may decrease by a factor of 2-3 as 
compared to the force without magnetic fields. At high 
Mach numbers, M. > A^^it, the drag force is slightly 
suppressed as compared to the unmagnetized case, but 
this reduction is more modest than for = 0°. Once 
again, at low Mach numbers {A4 < 0.75), the drag force 
is stronger than in the unmagnetized case. 

As already said, the component of the force perpendic- 
ular to the velocity of the perturber, fporp, will tend to 
induce a change in the direction of the velocity (note that 
wc force the body to move along a straight line). We will 
use the following sign convention for fporp- For an angle 
6 in the interval < O < 7r/2, i^pcrp > will mean that 
O > 0, in our convention. In Figure [T2| Fporp is shown 
as a function of Mach number. The magnitude of i^pcrp 
may be comparable to the drag force. For instance, at 
M = 1.2, the perpendicular force is only a factor of 2 
smaller for T = 1.41 and a factor of 3 for T = 1. Given a 
certain supersonic velocity, i^pcrp increases with T, while 
F]jF shows the opposite behaviour. For supersonic mo- 
tions with angle Q = 7r/4, i^porp is always positive and 
increases monotonically in time (see Fig. I13|) . This im- 
phes that Fperp will tend to redirect perturber 's velocity 
to a higher Q. For the cases shown in Figure [T^l the 
perpendicular force saturates in the run of the calcula- 
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tion (< — 40tcross) only in two cases; for Ai = 0.5 and 
T = 0.71, and for M = 0.75 and T = 1.41 (this latter 
case is shown in Fig. I13p . In some cases with subsonic 
Mach numbers, the perpendicular component of the force 
is initially positive, achieves a maximum and then stars 
a linear decline up to negative values (see Fig. [T^ . 

5.4. Dependence of the drag force on O 

In Figure [T3] we plot the drag force as a function of 8, 
for T = 1 and 1.41. The dependence of F^p on 8 is 
not always monotonic. The strongest variation of Fof 
with 8 occurs for A4 — 0.9. For this Mach number, the 
drag force may decrease by a factor of 2-3 from 8 = 
to 8 = 30°. For Ai = 1.2, the drag force may change up 
to a factor of 2 depending on the 8-value. For M = 0.3 
and 1.4, the drag force depends gently on the angle. 

In many astrophysical scenarios, the perturber will be 
subject to an external gravitational potential and will 
describe a nonrectilinear orbit. Sanchez-Salcedo & Bran- 
denburg (2001) numerically treated the orbital decay 
of a perturber in orbit around a unmagnetized gaseous 
sphere. They found that the "local approximation" , that 
is estimating the drag force at the present location of the 
perturber as if the medium were homogeneous but taking 
appropriately the Coulomb logarithm, is very successful. 
Consider now a perturber on a circular orbit in a magne- 
tized medium. If the orbit lies in a plane perpendicular to 
the magnetic field, the attack angle 8 is always 7r/2. In 
the local approximation, the maximum drag for Q — tt/2 
and T = occurs at 1 and at 1.7 for T = 1. 

However, if the plane of the circular orbit is parallel to 
the direction of the initial magnetic field, 8 will change 
periodically in time as 8 = 8o + ^ot. Therefore, if the 
local approximation is valid, one can estimate the mean 
drag force over a rotation period, which is approximately 
equivalent to take the mean value of Fof over 8. In par- 
ticular, for T = 1, the 8-average drag force is maximum 
at A4 « 1.4. This example illustrates how F^p may de- 
pend on T and on the inclination of the orbit respect to 
the magnetic field lines. A more detailed analysis of the 
drag force on a body on a circular orbit will be given 
somewhere else. 

6. SUMMARY AND DISCUSSION 

Understanding the nature of the DF force experienced 
by a gravitational object that moves against a mass den- 
sity background is of great importance to describe the 
evolution of gravitational systems. In this work, we in- 
vestigated the DF on a body moving in rectilinear tra- 
jectory through a gaseous medium with a magnetic field 
uniform on the scales considered. In linear theory, the 
problem is largely characterized by three dimensionless 
parameters, M, which is defined as the ratio of the par- 
ticle velocity to the sound speed of the uniform gas, T, 
defined as the ratio between the Alfven and sound speeds, 
and 8, the angle between the magnetic field direction and 
the particle velocity. We find that magnetic effects may 
alter the drag force, especially for T > 0.5, because the 
magnetic field affects the flow velocity field in the perpen- 
dicular direction of the ambient field lines, and thereby 
the morphology of the wake. Note that the plasma beta, 
defined as the ratio of gas to magnetic pressure, is 2/T^ 
for an isothermal system. 
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Fig. 14. — Dependence of the drag force on the angle for T = 1 
(upper panel) and T = 1.41 (lower panel). The numbers given at 
each curve represent the Mach number M. The drag force was 
computed at t = 40fcross • 

There are two major differences between the magne- 
tized and unmagnetized case. One conceptual difference 
is that, while gravitational focusing in a unmagnetized 
medium always generates a positive density enhance- 
ment, this is not the case in a magnetized medium (see, 
e.g.. Figs. [21 [3] and [9|). A second result is that the peak 
value of the drag force is not near A4 = 1 for a mass mov- 
ing in a magnetized medium. In fact, the sharp peak of 
Fup at Ai = 1 found in the T = case is no longer 
present in a magnetized medium with T > 0.5. For in- 
stance, for a perturber in perpendicular motion to the 
field lines (8 — tt/2) in a medium with T — 1.41, the 
drag force is essentially constant from A4 = 0.5 to 1.4 
and its maximum is located around Ai — 2 (see Fig. ITOl) . 
The flat plateau in the drag force between Ai = 0.5 and 
1.4 is partly because of the extra rigidity of the magnetic 
field in the x and y directions. 

For a body traveling along the field lines, i.e. 8 = 0, 
the steady-state problem can be treated analytically. We 
focus first on this case. For T 7^ 0, the drag force presents 
two local maxima (see Fig. [T]); one is located in the sub- 
sonic branch (at A^crit) and the other peak value is at 
the supersonic branch [at Ai — max(l, a/2T)]. When 
the velocity of the perturber is supersonic and super- 
Alfvenic (and 8 = 0), the DF force in a magnetized 
medium is weaker than it is in the unmagnetized case 
by a factor of (1 — rj), with 77 = (ca/Vq)^. The physi- 
cal reason is that the medium becomes more rigid in the 
radial direction and, hence, the opening aperture of the 
modified Mach cone is the same as that in a unmagne- 
tized medium with effective sound speed (c^ -I- c^)^/^, 
but the density enhancement is smaller by a factor of 
(1 — 77). By contrast, the drag force for subsonic veloc- 
ities is stronger if the medium is uniformly magnetized. 
For 8 = 0, an underdense region is formed upstream 
because of the gas channeling along the direction of the 
magnetic field, following the path of less resistance. The 
steady-state theory predicts that the gravitational drag 
on a body with 8 = vanishes at Mach numbers in 
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the following two ranges: (1) at Ai < Aiait and (2) 
at min(l,T) < A4 < max(l,T). However, using time- 
dependent analysis we find that the DF force asymptoti- 
cally approaches to a nonzero steady-state value at these 
Mach numbers. For T > 0.4 (still = 0), the drag force 
is maximized for perturbers moving at a Mach number 
close to TWcrit (Fig. [5]). At Mach numbers around A^crit, 
the density enhancement is large but negative in a cone in 
front of the body. At those Mach numbers, the DF may 
be even more efficient than in the stellar case. For exam- 
ple, for a medium with T — 1.41, the drag force peaks 
between = 0.6 and = 1.1. As a consequence of the 
stronger DF force, subsonic massive objects in a orbit 
elongated along the magnetic field lines in a constant- 
density core of a nonsingular gaseous sphere will suffer 
a orbital decay faster if the medium is pervased by a 
large-scale magnetic field. 

We have also explored the 9-dependence of the DF 
drag. For Mach numbers around A^crit, the drag force 
exhibits the strongest variations with & (see Fig. fTl)) . For 
magnetized media with T > 1 and regardless the exact 
value of 0, we find that (1) the drag force for subsonic 
perturbers is higher by a factor of 1.5-2 than it is for the 
unmagnetized DF drag, and (2) for supersonic perturbers 
{M > 1), the magnetized drag force is always weaker 
than the unmagnetized drag force. At intermediate Mach 
numbers, 1.1 < Ai < 1.4, the drag force is a factor of 2-3 
weaker than it is in the absence of magnetic fieldfl At 
high Mach numbers, M > {l.T + T^y/"^, the suppresion 
of the drag force is more important at small values of Q 
fFig. At these high Mach numbers and for an angle 
of Q — tt/2, the drag forces are similar with and without 
magnetic fields (Fig. [TU)) . 

As a consequence, supersonic massive objects may 



make their way more slowly to the center of the sys- 
tem if the medium is pervased by a large-scale mag- 
netic field. As a model problem, consider a singular 
isothermal spherical cloud threaded by a uniform mag- 
netic field and a small-scale random magnetic field with 
Alfven speed Ca everywhere constant. The density profile 
of the cloud is given by p{r) = (c^ +c^/2)/27rG'r^, where 
Cs is the isothermal sound speed. The circular speed 
is Vq = •\/2cf~fcf. Since the effective sound speed is 

•\/ 7C^ + 2Cq/3, the Mach number of a body on a quasi- 
circular orbit is 



M = 



2ci 



3^a 



1/2 



(52) 



which varies from 1.1 to 1.4 depending on the value of 
Ca and whether the perturbations are isothermal or adi- 
abaticQ. If the uniform magnetic field component has a 
T-value between 1 and 1.41, the time for the perturber's 
orbit to decay will be a factor of 2-5 larger than the 
corresponding decay time for T = 0. Our results demon- 
strate that, in the presence of ordered magnetic fields 
with T > 0.7, the role of the magnetic field on the drag 
force should be taken into account to have accurate esti- 
mates of the timescales of orbital decay via gravitational 
DF. 
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APPENDIX 

A. FOURIER TRANSFORMATION: AXISYMMETRIC CASE 
The three-dimensional Fourier transform of a perturbed variable ,f{r) is given by 



1 



In the Fourier space. Equations ([T5|) and p7)) are transformed into: 



(Al) 



(A2) 



(A3) 



where pp is the mass density of the perturber, thus V^^ — AirGpp. In order to have an equation for a, we will eliminate 
Pz- From Eq. (jA3|) . we have 



ci 



and substituting into Eq. (|A2p we find 



[M^kl - fc^l + (1 - M-^)T^)] --^(J^-^)Pr>^ 



(A4) 



(A5) 



^ This factor may be larger at later times because the magnetized 
drag force saturates, whereas it increases logarithmically in time 
in the unmagnetized case. See Figure |8] for an evolved stage. 

* In the nonmagnetic simulations of the orbital decay of a single 



black hole due to gaseous DF in Escala et al. (2004), the velocity of 
the black hole is initially supersonic (A4 = 1.4) and remains barely 
supersonic through most of the simulation. 



DYNAMICAL FRICTION 



19 



In the absence of magnetic fields (T = 0), the above equation reduces to 

iM^kl-e)a^~^Pp, (A6) 

and the standard steady-state equation for the wak e pa st a gravitating body is recovered. At velocities much larger 
than the Alfven speed, Ai ^ ca/cs — T, Equation (jA5|) is simplified to 



rPp- (A7) 



(i+T2)cr 

By comparing the above equation with Equation (|A6[) . we see that the response of the medium in this case is indis- 
tinguishable to that of an unmagnetized medium with sound speed (c^ -|- c\Y^'^ = 0^(1 + T^)^/^. 

It is interesting to note that when Vb = ca 7^ 0, the right-hand-side of Equation (|A5P vanishes and thereby the 
solution is a = 0, implying that the steady-state configuration satisfies V • t?' = 0. Obviously, the drag force is exactly 
zero in this configuration. 

There exist two situations where the differential equation (|A5I) is not well-posed: (1) dX M. = \ and (2) at the critical 
Mach number, A^crit, satisfying that 

1 + (1 - M-,l)T^ = 0. (A8) 

So that 

Merits (l + T-2)-'/^ (A9) 

It is clear that A^crit < 1. If the dynamics is dominated by the magnetic field, i.e. when T ^ 1, then A^crit — > 1. If 
not specified, we will consider M. ^1 and M ^ A^crit throughout this section. 

We will now calculate the solution of Eq. (IA5|) when the perturber is a point mass M, so that pp — Mj (27r)'^/^, which 
corresponds to the Fourier transformation of Pp(r) = M(5(r), to obtain the Green's function. Using the convolution 
theorem, it is possible to evaluate a for any general distribution p^. Hence, we solve for 



where 



(All) 



and 

e= 1 + (1-X-2)y2^ (^^2) 

The integral (|A10p along fc^ is evaluated by transforming to the complex plane. It is convenient to define 7^ = 1— A^^/^. 
Either if M. stands in the range < Al < AJcrit or in the range min (1, T) < AJ < max (1, T), then 7^ > and thus 
none of the poles lie on the real axis. Hence we may use Cauchy's residue theorem to obtain: 

a r = -^^ , A13 

For Mach numbers in any of the two ranges: Alcrit < A^ < min (1, T) and Al > max (1, T), the integrand has 
poles on the real axis. Hence we make 'indentations' in the contour at the position of the poles. We consider first 
Mach numbers larger than max(l,T). Then, for z > 0, we close the contour at -|-ioo, leaving the poles outside the 
contour to preserve causality, whereas for 2; < 0, we consider a domain containing the lower half-plane, that is where 
Im(A:2) < 0, and the contour slightly above the real axis, so that the two poles lie inside the contour. More specifically, 
for z < 0, the integration over can be evaluated as 

dk. = / , ..... rffc. = sin , (A14) 



A;2 - M\-^kl kl + kl- ^Ikl j,kn \ 71 

where 7J = — 7^ = — 1 and kj^ ~ k'^ + k^. The integration over k^ and ky can be carried out in polar coordinates 

(kx = kfi cos (j) and ky — fc^^sin^): 

'"sin(^^j e^'=«™^*rf,/. dkn 

/ sm Jo(fcfli?) dkR 

71 Jo V 71 / 

= {7^ (A15) 
otherwise. 
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At Mach numbers in the interval A^ciit < < niin(l, T), causality is perserved at 2; > if we consider a domain 
containing the upper half-plane and the two poles lie inside the contour. Putting all together, the solution is 



a(r) 



A(l - ri)GM 



1 



(A16) 



where 



2 if > max(l,T) and z/i? < -|7|; 
2 if TWcrit < < min(l,T) and z/R > \-f\; 
1 if < Xcrit or if min(l, T) < M < max(l, T); 
otherwise. 

This result can be compared with that in Dokuchaev (1964) by noting that he used the variable A to denote the 
combination (1 — ry)/^. Dokuchaev (1964) found the same functional form for a, but failed to divide correctly the cases 
according to the Mach number. In the absence of a background magnetic field, which corresponds to ^ = 1 and 77 = 0, 
we recover the classical form derived by previous authors. 

B. FOURIER TRANSFORMATION: PERTURBER'S VELOCITY PERPENDICULAR TO THE MAGNETIC FIELD 

We consider a gravitational perturber moving at constant velocity voy in an unperturbed medium with density po, 
sound speed Cg and a magnetic field Bo = BqZ. Therefore, the velocity of the perturber and the magnetic field are 
perpendicular. We are interested in the stationary wake formed behind the body. To do this, we assume that perturber 
is at rest and feels a wind with velocity —Voy at infinity. Following the same approach as in Appendix]^ the two 
coupled governing equations in this geometry are given by 



,& 



1 



1 

c? dz^ 



Tw- 



in the Fourier space, these equations are transformed into 



AttG 



dy"^ 



k^a + -irk^^ 
ci 



In order to have an equation for a, we eliminate using 



I3z = 



1 



AttG 



Pp + [M^kl - k'')a 



T2fc2 

and we obtain the solution of the density disturbance in Fourier space 

M 



(Bl) 
(B2) 

(B3) 
(B4) 

(B5) 



il + T')klk'-M'k^y-i-klk' 



Pp- 



(B6) 



Of course, the structure of a in this case is different than in the axisymmetric case (Eq. IA5|) . The inverse Fourier 
transform cannot be derived analytically but since the above equation is well-posed for any A^-value if T 7^ 0, we 
expect the drag force to be a continuous function of A4. 
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